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Abstract 

Collider data can play an important role in determining the parton distribution func- 
tions of the nucleon. I present a formalism which makes it possible to use next-to-leading 
order calculations in such a determination, while minimizing the amount of numerical 
computation required. 



Submitted to Nuclear Physics B 



^Laboratory of the Direction des Sciences de la Matiere of the Commissariat a I'Energie Atomique 
of France. 



1. Introduction 



Perturbative quantum chromodynamics gives an excellent description of short-distance scatter- 
ing processes at present-day colliders. The perturbative description relies on our ability to compute 
short-distance matrix elements in nonabelian gauge theories. It also relies on our understanding of 
factorization, which permits a separation of the process-dependent short-distance aspects from the 
universal, process-independent long-distance ones. The long-distance parts of scattering processes 
are captured in the parton distribution functions of the scattering nucleon(s). Along with the 
running coupling ccg, they are the only ingredients needed from outside perturbation theory for a 
description (up to subleading power corrections) of collider scattering processes. Precise knowledge 
of parton distributions is important in the quest for physics beyond the standard model. At current 
or planned high-energy colliders new physics must necessarily be detected against an omnipresent 
background of QCD or QCD-corrected events. 

The momentum evolution of parton distribution functions and of the running coupling is 
also governed by perturbative equations, so that it is only their values at a fixed scale which are 
required inputs from outside perturbation theory. Such input parton distributions may someday be 
calculated on the lattice or by other nonperturbative means, but at present they must be extracted 
from experiments. The modern approaches [1,2] involve global fits of ncxt-to-leading order theory 
to all available experiments. The experiments involve different scale arguments to the distribution 
functions, but as these arc related by the above-mentioned perturbative evolution equation, we can 
regard the fits as determining the distributions at a certain fixed scale Qq. 

To date, it is primarily deeply-inelastic scattering data that has been used in the global fits. 
(There is some ad-hoc use of collider data; for example, MRS [1] make use of certain points from 
the lepton asymmetry distribution from CDF and DO, relying on the fact that NLO corrections 
are small for these points, and on the availability of an analytic calculation for this quantity.) Yet 
there is a wealth of collider data which may yield important constraints on the parton densities, in 
particular the gluon distribution. The latter enters into DIS calculations only at higher order in 
the coupling, but is a dominant contribution in hadron-hadron collisions. 

The deeply-inelastic structure functions used in the fits can be described by a convolution of an 
analytically-known function with the parton distributions. This makes their use in a fit computa- 
tionally feasible. In contrast, for next-to-leading order differential cross sections at colliders (both 
hadron-hadron and lepton-hadron) , even the short-distance part must necessarily be calculated 
numerically. This is a consequence of the relatively complicated structure of phase space, once one 
allows for arbitrary experimental cuts and jet algorithm. 

In principle, existing NLO jet programs can be used in the global fits. Were we to try to use 
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them, however, we would have to regenerate the parton distributions (for example, in the form of 
numerical tables) at each iteration of a fit, and then re-run the jet program given the new values. 
This latter part of this procedure would be extremely time-consuming, and completely impractical 
for fits requiring more than a handful of iterations. 

Yet much of the calculation of a jet differential cross section — the jet algorithm and experimen- 
tal cuts, the numerical integration over real emission, the balancing of real and virtual contributions 
— is in fact independent of the precise form of the parton distribution functions. Thus it would be 
useful to re-organize the calculation so as to minimize the amount of computational work needed 
at each iteration of a fit to experimental data. 

One might imagine pre-generating values for a grid of parameter values, and interpolating 
between them to find the best fit, but if in a realistic case one has 15 parameters, and wants (say) 
to consider ranges discretized into 10 values each, one would have to generate 10^^ different points, 
clearly a hopeless task. 

Fortunately, there is a better way; read on to discover it. Graudenz ct al. [3] have previ- 
ously presented an approach to using certain jet distributions from lepton-proton scattering. The 
approach presented here has certain elements in common with theirs (such as the use of Mellin 
transforms), but is fully general. It can be applied to arbitrary differential distributions in both 
lepton-hadron and hadron-hadron scattering. It is also free of certain theoretical restrictions and 
numerical limitations present in their approach. 



2. A Prototype 

Let us first consider leading-order calculations for a glueball G in a quarkless version of QCD. 
This will serve as a warm-up exercise for the real-world case. 

In this case we have only one distribution to consider, fg^G{x,Q^)- The total n-jet cross 
section, subject to experimental cuts, in glueball-glueball scattering, is given by 

o'n = / / dxidx2 / dLlPS{xikG + X2kQ ^ {kiji^-^) 
Jo Jo J 

X fg^G{xi,I^F{{ki},Xl,2))fg^G{x2,IJ'%{{ki},Xi^2)) (2-1) 

X a"(/^l({fci}>a;i,2))o-(55 {ki}) Jn^n{{ki}) , 

where LIPS stands for the Lorentz-invariant phase-space measure. In this equation, a stands for 
the usual leading-order partonic differential cross section with the running coupling ag set to 1, and 
fg^oix, M^) is the gluon distribution inside the glueball. Note that the ki are implicitly dependent 



on xi and X2 as well. The renormalization and factorization scales fiR and fxp — typically something 
like a jet Et — also depend on xi^2 and the final-state momenta (thereby violating, for example, 
one of the theoretical restrictions in the work of Graudenz et al. [3]). 

The jet algorithm is represented by Jn<-„, which evaluates to 1 if the original n-parton con- 
figuration yields n jets satisfying the experimental cuts, and otherwise. The precise form of the 
jet algorithm is not important for the formalism presented here. 

The parton distribution functions satisfy evolution equations, whose solutions can be written 
using a universal evolution operator E{x,as{Q^).,ao). (The initial coupling ao = (^s{Qo) one of 
the parameters we will want to fit.) Explicit forms for the evolution operator can be found in the 
paper of Furmanski and Petronzio [4] and elsewhere. The solutions can be written in the form 

fg^cix, Q') = E{x, a,{Q^), ao) ® fg^cix, Ql) , (2.2) 



where 



A{x) (g) B{x) = [A0B] (x) = [ dy [ dz S{x - yz)A{y)B{z) (2.3) 

Jo Jo 



defines the convolution symbol (g). 

A Mellin transformation turns these convolutions into multiplications, 

/I^g(Q') = E%a,iQ'),ao)f^^GiQl) , (2-4) 

in which 



/' 

Jo 



dx x^-M(x) (2.5) 



is the Mellin transform of A. 

To recreate the original function, use the inverse Mellin transform. 



f C+200 



A{x) = ^ f dz x-^A^ . (2.6) 

The contour should be chosen to the right of all the singularities of A^ . 

We can thus write 

1 j'C+ioo 

fg^G{x,Q') = — / dz x-^E^{a,{Q%ao)fl^Gi.Ql) ■ (2-7) 

All the parameters (except ao) that we wish to fit are contained in fg^Q{Q'^). This latter function 
is independent of all integration variables except z, and thus can be pulled out of the numerical 
integrations in eqn. (2.1). The remaining zi^-z contour integrals are to be performed during the 
fitting procedure, but this involves only a double sum of the gluon distribution function multiplied 
by precomputed numerical coefficients. We will be able to make use of techniques described in 



ref. [5] to find numerically efficient contours for the Zi integrations. I postpone the discussion of 
the choice of contours, and the method for performing the contour integrals, to a later section. 
Thus at each step of the fitting procedure, we must compute only 

1 nc+ica nc+ioo 

-J^ dz^fluiQDflLciQD^''''^ (2-8) 

J c—ioo J c— zoo 

where S^i'^^ g^j-g precomputed coefficients given by 

^zuz^^l f ^^^^^2 [ dLlPS{xikG + X2k'G ^ {ki}2=i) 
Jo Jo J 

X x^^'x^'^E'^ {a,{fj,%{{ki}, xi,2)), ao)E'' ia,{p,l{{h}, ^1,2)), ao) ^^"^^ 

X Ol'iil4ii{ki},Xi^2))&{gg {ki}) Jn^ni{ki}) ■ 

Since S contains all of the short-distance process-specific dynamics, but is independent of the 
parton distributions (and indeed, of the nature of the parent hadron), it is appropriate to call it a 
universal cross section, in this case to leading order in quarkless QCD. 

This procedure works for the parameters in fg^Q^Q^) because T, is independent of them. It 
does not work for oq, the remaining fit parameter, because cross sections and distributions depend 
on Q!o in a (complicated) non-polynomial fashion. There are several approaches we can take here. 
The simplest is to generate the coefficients E^^'^^ for a set of ao around a "canonical" value (e.g. 
as{M^) = 0.117), and then fit using interpolation. While this approach would be vastly too 
time-consuming for a large number of parameters, it is acceptable for a lone parameter. 

In general, we don't want to fit total cross sections, but rather differential distributions. The 
above discussions go through just as well for the latter. 

3. Leading-Order Fits 

Let us turn next to leading-order fits in the real world, specializing to the Tevatron. The 
formulae are quite similar to those in the previous section, except that we must sprinkle a variety 
of indices in appropriate places. 

The ra-jet differential cross section (in the variable X) is now given by 



/ / dxidx2 / dLIPS(xikp + X2kp ^ {kij^^i) 
„h Jo Jo J 

X fa^p{xi,nU{h}, xi,2))fb^p{x2,ti%{{ki}, xi^2)) S{X - X {{ki})) (3-1) 



dX 

ao 



X <{f^li{ki},xi,2))a'^''{ab^{ki}) Jn^^Uki}) . 
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The partonic cross section u is implicitly summed over all different possible final states, and the 
partonic types a, h are summed over the gluon and all relevant quark and antiquark distributions. 
Eqn. (2.4) is replaced by 



jl^^{Q^) = El,{a,{Q\ a^)n^^m ■ (3-2) 

The matrix is most easily expressed in a basis of the evolution eigcndistributions q — q, q + q, 
along with gluon distribution g and the quark singlet distribution S. The initial parton distributions 
are often expressed in this basis as well, but we need the evolved distributions in the usual flavor 
basis. It is therefore convenient to use in a form where its left index is in the flavor basis, and 
its right index is in the eigendistribution basis. 

Each step of the fitting procedure involves the computation of 



J c — ioo *J c — ioG 



(with implicit summation over a, 6), where 

— = / / '^^i'^^^ dLIPS{xikp + X2kp ^ {ki}1^^ 

xx-"^a;2-'^S:,^„(a.(MF({A;a, 2:1,2)), «o)^fe4K(^l({^J' ^1,2)), ao) ^^'^^ 
X <(/x2j({A;,},xi,2))aL°(a'6' ^ {h}) Jn^n{{ki})S{X - X{{ki})) . 

In a numerical computation, the delta function is implemented by binning in X. 



4. Next-to-Leading Order Fits 

At ncxt-to-lcading order, the structure of the cross section is more complicated; we must com- 
bine virtual corrections with real-emission ones. This is a delicate procedure, because each of these 
contributions is independently infrared divergent, and only their sum is well-defined. Moreover, 
from a practical point of view, the only practical infrared regulator is dimensional regularization, 
which does not mesh naturally with numerical calculations. The notion of a parton resolution — 
most simply an invariant mass Smim with a pair of partons unresolvable if (fcj -|- kj)'^ < Smin 
— offers a safe passage through these treacherous divergences and cancellations. The basic idea 
is to evaluate analytically the infrared-divergent contributions from real emission over the phase 
space for (color-adjacent) unresolved partons, and to add this contribution to the virtual correc- 
tions. The poles in the dimensional regulator e then vanish, and one can take the four-dimensional 
limit of all expressions. The remaining contributions, over the phase-space for resolved partons, 
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are computed numerically. Upon adding the two contributions, the dependence on the resolution 
parameter disappears (in the limit Smin — ^ 0) and one obtains the next-to-leading order differential 
cross section. Color ordering plays an important role because the notion of resolved or unresolved 
partons is defined independently for each color permutation. 

Several variants of this approach to general processes in collider physics have been presented 
in the literature, in particular the 'slicing' [6,7] and 'subtraction' [8,9] ones. The details are again 
not important for our purposes. The exposition below assumes the use of either the pure slicing 
or a 'restricted subtraction' method, for which a subtraction of singular pieces is performed only 
inside the region Sging < Smin- Analogous formulae can however be written down for the subtraction 
method of refs. [8,9]. I shall indicate the restriction to resolved configurations via an additional 
superscript R on, and an additional argument Smin to, the appropriate leading-order partonic cross- 
section. Similarly, a superscript U will denote the restriction to unresolved configurations. 

Thus 

dX dX dX dX ' 

The contribution of configurations with n + 1 resolved partons is. 



NLO,R /•! /•! 

dX 



dSa^ 



V / / dxidx2 f dLIPS{xikp + X2kp ^ {ki}'l+^) 

X /a^p(a;i,/x|({^i},xi,2))/6^p(x2,/x|({/ci},xi,2))(5(X (4-2) 

X <+H/x|j({/Ci},Xi,2))(T^°'^(a6^ {A:J;Sniin) Jn^n+l{{ki}) . 

In this equation, Jn^„_|_i({A;i}) evaluates to one if a given (n + l)-parton configuration yields n 
detected jets, and vanishes otherwise. It is crucial that both the jet algorithm and the differential 
cross section under consideration be infrared-safe, to wit the treatment of an event must not change 
under the addition of an arbitrarily soft gluon, or the splitting of any parton into two collinear 
partons. Otherwise, the details of the jet algorithm are again unimportant for the formalism 
presented here. 

What is more important for us is the structure of the unresolved pieces, 

dX dX dX ' ^ ' ' 



The first term combines the virtual corrections with the singular integrals over the unresolved 
phase space, in a crossing-invariant fashion. This term would be present whether or not the initial 
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dX 



state contained colored partons, 

ab 

X fa^p{xi, l4p{{ki}, Xi^2))fb^p{x2, fljpiih}, Xl,2)) S{X - X{{ki})) 



= / / dxidx2 [ dLIPS{xikp + X2kp ^ {ki}i^-i^) 
Jo Jo J 



(4.4) 



The second term crosses colored partons from the final state to the initial state; it is thus 
absent in e+e~ scattering. In our case, it takes the form 

NLO,C A /•! 



d6a 



I / dxidx2 dLIPS{xikp + X2kp ^ {ki}'^^-^) 



X [Ca^p{xi, H%{{ki}, Xi^2y, Sn,in)fb^p{x2, f^piiki}, Xi,2)) 

+ fa^p{xi,IJ,%{{ki}, Xl^2))Cb^p{x2, IJ,%{{ki}, 2:1,2); Smin)] S{X - X{{ki})) 

X 2^1,2))^^° (a6 ^ {ki}) Jn^n{{k^}) , 

(4.5) 

where Ca^p is a crossing function as introduced in rcf. [7]. The crossing functions are factorization- 
scheme dependent; they can be expressed in terms of scheme- independent functions A^^p and 
scheme-dependent functions Ba^p as follows, 



with 



Aa^p{x, Q2) = Kf^j,(x) fb^pix, Q2) , 



(4.6) 



(4.7) 



Ba^p{x, Q2) = Kf^b(x) ® h^p{x, Q^) . 

Expressions for the kernels K^'^ are given in ref. [7], and expressions for their Mellin moments 
K^'^ and K^'^ in ref. [5]. 

With these moments, we can define 



K 



C,z 
a*—b 



N 
2^ 



T^A,z 1 / •''mill 



I ^a^b 



this allows us to write 



and then using eqn. (3.2), 



ci^p{Q')=K':Cbfb^pm, 



(4.8) 



(4.9) 



(4.10) 



Each step of the fitting procedure involves the computation of the same quantity as in eqn. (3.3), 
where now 

dX ~dX ^ 'dX ^ 'dX 

~ IX '^'dX '^IX ^~dX 

The first term, dTi^^'^^'^'^ / dX , is given by eqn. (3.4); the latter quantities are given by the 
following equations, 

= dx,dX2 J dLlPS{x,kp + X2kp^{ki}^+^')6{X -X{{ki})) 

X x-'^X2''EXMf4i{ki}, xi,2)), ao)E',?t,ia,i,xl{{h}, ^1,2)), ao) ^"^'^^^ 
X <+'(/x|({A;i},xi,2))aL°'^(a'6' ^ {ki}-s^i^) Jn^n+ii{ki}) , 



^^^NLO,F:2.,22 ,1 ,1 

"IF 



/ / da;i(ia;2 / dLIPS{xikp + X2kp ^ {A^^j^^i) (5(X - ^({A:^})) 
7o Jo J 

X ar-^^ar2 ^^£;^?Ja«(/x|({A;i}, ^1,2)), ao)E^?^MfiU{ki}, ^1,2)), ao) ^^'-^^^ 
X <+'(/x|({fci},xi,2))<5a^^°(a'6' ^ {A;i};s^i„) J„^„({A;i}) , 



^NLO,C:z„22 ,1 ,1 



/ / dxidx2 f dLIPS{xikp + X2kp ^ {ki}^^^) 6{X - X{{ki})) 
Jo Jo J 

K^C\KlMtiU{ki},xi,2)),ao)E^,?,iaMa^^^ 



V — 21™-Z2 



+EXiasipli{k},x,,2)),ao)K^/^lEl^^{as{^^^^^^ 

X <+HMli({Ai},Xi,2))aL°(a'6' ^ {fci}) Jn^niih}) . 

(4.14) 

The calculation of dS^^'^^/dX for a given (^1,^2) amounts to doing the usual next-to-leading 
order calculation, with the structure functions fa'^p{x,Q'^) replaced by x~^E^,^ (a remaining a 
free index), and the crossing functions Ca'^p{x,Q^; Sram) replaced by x~^K^,'^g^Ef^. 



5. Integration Contours 



The contour integrals in eqn. (3.3) must be performed numerically. This involves several 
distinct choices. First, we must choose an integration contour; then, we must choose a set of points 
along the integration contour, or in the case of hadron-hadron scattering, a set of points in the 
plane defined by the product of the two contours. These points must be chosen in advance, of 
course, so that we can evaluate the S matrices at them. How should we make these choices? 
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The simplest choice is to pick the textbook contour, parallel to the imaginary axis, but displaced 
to the right. This is what was done in the old papers of Gluck and Reya [10] on parton evolution, 
and also in the paper by Graudenz et al. [3]. The integral evaluated was actually j^-iT ^ relying 
on the fall-off of the integrand as T — ^ oo to drop the remaining terms. This finite-interval integral 
was evaluated using Gaussian quadrature. As discussed by Gluck, Reya, and Vogt [11] and in 
refs. [3,5], however, the integral falls off rather slowly in this direction, only as a power of the 
integration variable z. For parton evolution, the pole structure is such that one can freely deform 
the contour into the left-half plane, whereupon the integrand falls off exponentially, improving the 
convergence of a numerical evaluation. Graudenz et al. [3] used the Mellin moments of the short- 
distance cross section in their formalism; these do not fall off in the left-half plane, and thus they 
were not able to deform the contour. (In contrast, the numerical Mellin inversions used nowadays 
in parton evolution programs do utilize a rotated contour [11].) 

In the formalism presented in previous sections, however, the Mellin moments of the short- 
distance cross section appear nowhere. We may note that the experimental cuts (on rapidity and jet 
transverse energy) effectively impose a minimum on the parton momentum fractions X\^2^ so that 
dS^i'^^/dX has no poles in the right-hand plane. Furthermore, since < xi^2 < 1; this quantity 
will also fall off exponentially as Zi oo, so long as Re 2, < 0. Thus just as in refs. [11,5], we can 
freely shift the contour into the left-half plane, so long as we stay away from the poles of fa^piQo) 
along the real axis. Indeed, much of the formalism developed in ref. [5] carries over to the choice 
of contours in the present paper. 

The desired contour would be the contour of steepest descent (which because of analyticity is 
also the contour of stationary phase, upon which our integrand will be purely real). For performing 
the inverse Mellin transform integrals required for the evolution of parton distributions, ref. [5] 
shows that it is possible to construct simple but good analytic approximations to such contours, 
and we may hope that those findings carry over to the contours required for evaluating eqn. (3.3). 

In any event, we should expect the contour of steepest descent for the leading-order calculation 
to be very close to that for a ncxt-to-lcading-ordcr calculation; furthermore, wc should expect it to 
not be very sensitive to the form of or parameters in the initial parton distributions. In practice, 
we have a reasonable idea of the initial distributions (the starting point for the fit, say an existing 
distribution set, will not be radically different from the best fit), so a good approximation to the 
contour of steepest descent we seek will be given by the contour of steepest descent in eqn. (3.3) 
for a leading-order calculation with any existing pdf set. 

It will be helpful to examine first the case of lepton-hadron scattering, which is simpler because 
there is only one inverse Mellin transform to compute (and hence a contour in only one variable 
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to choose). I shall discuss the technically more complicated hadron-hadron case in the following 
section. We can set up the fitting procedure following the same formalism developed in previous 
sections. At leading order, each iteration of a fit to a n-jet differential cross section in deeply 
inelastic scattering (the count excludes the remnant) requires computing 



1 r+^°° . ,^2,d^i 



dzf^^M):iv'' (5-1) 



(with implicit summation over a), where 



— « = y dxj dLlFSixkp + k, ^ K,{k,}U) x-^E^^,^iasifili{h},x)),ao) 



(5.2) 



In principle, we could let our chosen contour depend on the variable X, but this is both 
cumbersome and computationally more expensive. We would rather find a single contour for all 
values of X, so long as we find a contour (along with a set of point along that contour) that yields 
an accurate evaluation of the integral for all X where we have data with small statistical and 
systematic errors. We could just take the cross section (the integral over X) but if X depends on 
the transverse energy Et — that is, if the distribution is differential with respect to the transverse 
energy as well as other variables — then the total cross section will be dominated by events just 
above the lower Et boundary. This happens because the cross section will be falling rapidly as 
a function of Et- It is probably better, for purposes of determining the contour, to weight the 
integral over X by a 'stratifying' function S{X) which compensates for the rapid decrease as a 
function of transverse energy. Take 



dXS{X)^\ (5.3) 



where fal^p{Qo) denotes the initial parton distribution at the start of the fit procedure. We 
can now follow the approach outlined in ref. [5], the only change being that derivatives of F 
must be computed numerically rather than analytically. We can, however, write down closed-form 
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expressions for the required derivatives, 



F'{z) = j dX S{X) j^dx j dUPS{xkp + k'^, {itj^^i) 



X X 



- E^,^{as{lJ,%{{ki},x)),ao) + — £;^,„(a<,(//|({fei}, x)),ao) 



F"{z) = J dX Six) dx J dUPS{xkp + fee ^ Kdki}i=i) 



X X 



In xE^,^{asin%i{ki},x)),ao) - 2\nx-^E^,^{as{n%{{ki}, x)), ao) 



dz 



X ar\lil{k',, {h}, xi,2))aL°(ea' ^ fc^, {h}) Jn^niih}) 5{X - X{k',, {h})) , 
dX S{X) dx J dUFSixkp + k^ ^ k'^, {/ci}"=i) 



XX - -ln^a;F^,„(as(/x|({A;j},a;)),Q;o) + 31n^a;^£^^,„(as(/x|({A;i},a;)),Q!o) 

-^^^x—E^,^{as{li%{{ki},x)),ao) + ■^E^,^{as{iJ,%{{ki},x)),ao) 

X ar\t^lik',,{h},xi,2))a'^''iea' ^ k'„{h}) Jn^ni{h}) SiX - Xiki,{h})) . 

(5.4) 

The procedure there can be summarized as follows. First, find the minimum cq of F along the 
real axis. Next define the curve 



z{u) = Co + ic2\/u + clc3u/2 , 



(5.5) 



where 



2F(co) 
V F"{co) ' 

3F"(co) ■ 

To evaluate eqn. (5.1), evaluate the transformed version. 



C2 



C3 



C2_ 

2tt 



f 

Jo 



du 



Re 



u 



e" (1 - ^C2C3V^) /,^t"^(Q2 



using a generalized Gauss-Laguerre quadrature. 



(5.6) 



(5.7) 



(5.8) 
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where the u'j are the zeros of the generahzed Laguerre polynomial Ln ^^^\u), and the weights are 
given by standard formulae [12], 



It;,- 



r(n+l/2) 



n! (n + 1) 



r(-i/2) 



K) 



(5.9) 



The number of points n required for evaluating eqn. (5.7) to a given desired precision, using the 
Gauss-Laguerre quadrature formula (5.8), remains to be investigated numerically, but experience 
with parton evolution [5] suggests that n ~ 5 should be sufficient for < 1% error. 



6. Contours in Hadron-Hadron Collisions 



In hadron-hadron collisions, at each iteration of a fit, we have not a single contour integral to 
perform, but the double integral (3.3). (I again use the LO cross section to determine the contours.) 
Defining 

dXS{X)—^^ ' ^^-^^ 

we have to find not an approximation to a contour of steepest descent, but an approximation to a 
'surface of steepest descent'. We can take this surface to be invariant under conjugation of each Zi 
separately, and use this symmetry to rewrite the contour-determining integral, 

-j^ nc+i(x> /-c+ioo 

— -r^ I / [dZi -F'(2l,'22) - (^^l <^^2 -^(^1,^2) - (i^l <i22 -^(^1, -^2) + <^22 -^"(^1, ^2)] 

47r2 7o Jo 

/-C+ioo /-C+ioo 

= -^/ / 'R.e[dzidz2 F{zi,Z2) - dzidz2 F{zi,Z2)] 

(6.2) 

It is not necessarily symmetric in zi ^ Z2, because the beams, detector, or observable (such as a 
parity-violating one) may not satisfy the required xi ^ X2 symmetry. 

The poles in the Zi lie along the negative real axis, so wc can freely deform the contours into 
the left-half plane. Now just follow the procedure used in the single-contour case. First find the 
minimum of F{zi,Z2) for real 2:1^2, and label its coordinates (ci,C2). Parametrize the surface using 
two variables ti^2, 

zi = ci+ iti + aiit\ + 012*1*2 , 



^2 = C2 + it2 + 022*2 + «21*1*2 , 
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(6.3) 



Expand the equation Im.F{zi{ti,t2), Z2{ti,t2)) = to obtain 

i?(3,0) 



On 



6F(2'0) 

ir(0,3) 



ai2 = - 



021 = - 



1 

6F(o.2) + ^ 

i?(3,0) ;^ 

+ 



6F(2.o) 6J 

i?(0,3) 



^(0,3) ^(2,0) _ 3^(1,2) ^(1,1) _^ 3^(2,1) ^(0,2) _ 
^(3,0) ^(0,2) _ 3^(2,1) ^(1,1) _,. 3^(1,2) ^(2,0) _ 



p{3,0) p{0,2) 

_p(0,3)^(2,0)_p(l,l) 



022 



6F(0'2) 

where I use the abbreviated notation 

_p0'l,j2) _ 



J = det 



i?(2,0) ^(1,1) 
^(1,1) ^(0,2) 



(2l,22) = (ci,C2) 



We expect the function to go hke F(ci, C2)e 9{ti,t2) ^ with 



i?(2,o) 



tlt2 + 



i?(2.0) 



2F(ci,C2)''' ■ F(ci,C2)-'^^ ' 2F(C1,C2) ' 
which suggests the change of variables 



'F(c.,c.)(A + A) Ffa...)(A-A) 



A_A 



t2 



'F(ci,C2)(A-A) F(ci,C2)(A + A) 

-Viii + 1/ ^ — 1 V'"2 



A+A 



A_A 



where 



A± = 



^(2,0) ^ ^(0,2) ^ ^(J7(2,0) _ i7(0,2))2 ^ ^ 



are the eigenvalues of the Hessian matrix at (ci, C2), and 



A = A+ - A_ 



(6.4) 



(6.5) 



(6.6) 



(6.7) 



(6.8) 
(6.9) 

(6.10) 



Using the changes of variables (6.3,6.7), we can rewrite the integral (3.3) to be performed at 
each iteration of a fit as follows. 



00 /"OO 



1 F(C1,C2) 

27r2 Jo Jo 



duidu2 



-/uIU2 



h{ui,U2) , 



(6.11) 
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with 



h{ui,U2) = Re 



(l - z(2aii + a2i)ti(ui) - i{ai2 + 2a22)t2{ui) 

-2aiia2itl{ui) - iaiia22ti{ui)t2{ui) - 2012022*2 (^i)) (6-12) 



Zl{ui),Z2{ui) 



{zi^2 and ti^2 are functions of ni^2 via eqns. (6.3,6.7)) 

We will again use the generalized Gauss-Laguerre quadrature formula, here in each of the 
variables Uj. This would ordinarily lead to using points having chosen the n-point formula; 
we may note, however, that points at the far end of the square away from the origin will give a 
negligible contribution, so that we can restrict the sum to points {u^,u^) with + < + u^, 



duidu[ 



^/UlU2 



2 g — Ul— U2 



ii>j2=i 



^ Wj,Wj^h{u°^,u%). 



31 '32 — ^ 
«0 +«0 <uO+v.O 



(6.13) 

(The n° are again the zeros of the generalized Laguerre polynomial Ln ^^^^ (u) , and the weights Wj 
are given by eqn. (5.9).) 

The number of points needed remains to be investigated numerically. Formulae for the various 
derivatives used above can be written for computer evaluation in the same fashion as in the case 
of deeply inelastic scattering, eqn. (5.4), 



di^zidi^Z2 



zi,Z2) = fat'p{Ql)flCl{Ql) J dX S{X) j'^ dxrdx2 j dLIPS(xiA;p + X2kp ^ {fci}r=i) 



d^^ 



dz^^ 
d^' 

X 



dz^^ 



^K'ai(^A^J'Fi{k^},Xl,2)),ao)] 

b~''-^a'a(«s(Ml({fci},a;i,2)),ao)] 



Z=Z2 



X <(Ai^({A;i},xi,2))a^"(a'6' ^ {h}) Jn^^{{h}) 5{X - X{{ki})) 

(6.14) 

where fal^p{Qo) denotes (as in section 5) the initial parton distribution at the start of the fit 
procedure. 

In a practical application of the formalism presented here, one would proceed as follows. One 
would first determine the contours appropriate to each distribution one wanted to use in (say) 
fitting the parton distribution functions of the proton. (It is plausible that different distributions 
could make use of a common contour, but this is not guaranteed.) One would then determine the 
required number of points in the inversion procedure. These computations would be done using 
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leading-order matrix elements, as dicussed above. Having chosen the points Zi^2, one would then 
perform a next-to-leading order computation of the quantities £;5]^o.nlo,f,nlo,c,nlo,r^^^ defined 
in section 4. One typically uses functional forms for the initial parton distributions whose Mellin 
transforms are known analytically as a function of the parameters and of the Mellin variable z. 

One would seek to minimize the of the fit to experimental data by an iterative procedure. 
The Mellin transforms of the initial parton densities may be evaluated numerically for a given choice 
of the parameters. The quantity h defined in eqn. (6.12) is then just a numerical inner product, 
with the 'metric' given by the dT:ab/dX matrices evaluated previously. At each iteration of such a 
procedure, one would use eqn. (6.13) to compute the desired observable. Only the numerical values 
of the Mellin transforms of the initial parton densities would need to be evaluated anew at each 
iteration of the fit. 



7. Conclusions 



Jet data collected at the Tevatron and at HERA can play an important role in determining 
the parton distribution functions of the nucleon. To minimize renormalization-scale uncertainties, 
and to take full advantage of the data, they should be fitted using next-to-leading order (or higher 
order) theory. The formalism presented in this paper makes it computationally practical to fit to 
distributions culled from jet data, using a modern next-to-leading calculation. It reorganizes the 
calculation so that most of the time-consuming computations in an NLO program are done only 
once, and that each iteration of a fit involves only the recomputation of the Mellin transform of the 
initial parton densities at a few points, and a weighted sum over those points. The number of points 
required can also be minimized via choice of a quadratic contour, using the same approach detailed 
in ref. [5]. A similar approach could be used to reorganize the calculation of next-to- leading order 
corrections to allow the extraction of parton-to-hadron fragmentation functions from collider data. 
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